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Abstract 

We study the evolution of multiple supernova (SN) explosions inside a pre-exiting 
cavity blown by winds from massive progenitor stars. Hydrodynamic simulations in 
one-dimensional spherical geometry, including radiative cooling and thermal con- 
duction, are carried out to follow first the development of the wind-blown bubble 
during the main sequence and then the evolution of the SN-driven bubble. We find 
the size and mass of the SN-driven bubble shell depend on the structure of the 
pre-existing wind bubble as well as the SN explosion energy E$n (= A r 5Arl0 51 ergs). 
The hot cavity inside the bubble is 2-3 times bigger in volume and hotter than 
that of a bubble created by SNe exploded in a uniform interstellar medium (ISM). 
For an association with 10 massive stars in the average ISM, the SN-driven shell 
has an outer radius of R ss « (85pc)-/V<^ and a mass of M ss (10 4,8 M Q )N^ at 
10 6 years after the explosion. By that time most of the explosion energy is lost via 
radiative cooling, while <^ 10% remains as kinetic energy and ~10% as thermal en- 
ergy. We also calculate the total integrated spectrum of diffuse radiation emitted by 
the shock-heated gas of the SN bubble. Total number of H Lymann-limit photons 
scales roughly as $13.6 ~ 10 61 Nsn and those photons carry away 20 - 55 % of the 
explosion energy. For the models with 0.1 solar metalicity, the radiative energy loss 
is smaller and the fraction of non-ionizing photons is larger, compared to those with 
solar metalicity. We conclude the photoionization /heating by diffuse radiation is the 
most dominant form of feedback from SN explosions into the surrounding medium. 
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1 Introduction 



Massive stars inject a significant amount of energy into the surrounding medium 
in various forms: photoionization/heating by stellar radiation and mechanical 
energy through stellar winds and supernova (SN) explosions, some of which 
later transforms into gas thermal energy and radiation. For example, OB stars 
deposit stellar radiation energy into the interstellar medium (ISM) with the 
luminosity of 10 4 < L^/L Q < 10 6 and the H Lymann - limit photon luminosity 
(hu > 13.6eV) of KFs" 1 < Q H < KPs" 1 flSmithl . l2006h . So the radiative 
feedback of an early type star during the main sequence (tms ~ 10 6 yrs) may 
amount up to ~ 10 51_52 ergs, which is comparable to typical SN explosion 
energy. 



These stars also inject mechanical energy into the ISM via stellar winds during 
the main sequence (MS) and the red supergiant (RSG) or luminous blue vari- 
able (LBV) phase, and also Wolf-Ra yet stage in the case of very massive stars 



with M* > 25 M Q 
of D„ : 



(jChu et all 12004 ) . The MS winds have a typical wind speed 

1000-2500 km s" 1 and a mass loss rate of M w ~ 10~ 7 — 10~ 6 M^ — _1 

i 



1© V r ~ 



mass loss rate, M ~ 10" 



10" 4 M 



yr~ 



(Chu et al.. 


2004; 


Dwarkadas. 


2005 



Smith! . 120061 ) . Owing to much longer duration the MS wind generates a bigger 
bubble, while the RSG wind produces a much smaller but denser shell inside 
the MS bubble created earlier. The m echanical ene rgy luminosity of the MS 
winds ranges 100 < L wind /L Q < 3000 flSmithl . l2006h . 



Massive stars with mass M > 8M end their lives as core-collapse SNe and 
deposit typically E$n ~ 10 51 ergs in the form of kinetic energy of the ejecta. 
Initially the SN ejecta expands freely into the wind-blown cavity and later 
interacts with the shell produced previously by the winds, generating various 
shocks, e.g., transmitted and reflected shocks. Most of the explosion energy 
is thermalized via such shocks. After the shock-heated gas cools adiabatically 
down to 10 6 K or after the forward ejecta shock catches up with the wind 
bubble shell, the supernova remnant (SNR) begins to lose a significant fraction 
of thermal energy via radiative cooling. The final energy budget deposited 
into the surrounding medium in the forms of thermal, kinetic, and radiative 
energy depends on the detailed history of radiative cooling and heat transport 
by thermal conduction, which in turn depends on the size and structure of the 
preexisting wind bubble as well as on the ambient ISM. For example, if a SN 
explodes in a dense environment, the explosion energy would be lost mostly by 
radiative cool i ng and the impact to the ISM should be minimal. According to 
Tenorio et al.l (jl990bl ). in the case of SN explosions inside wind-blown bubbles 
in the typical ISM with nn,o = 1 cm -3 , the remaining kinetic energy fraction 
is Ek ~ 0.3Esn at the end of SNR evolution, while the remaining thermal 
energy fraction, E th ~ 0.2 — 0AE SN . 
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Feedback from massive stars should play a significant role in the formation and 
evolution of galaxies. Especially, feedback from the first-generation SNe affects 
the formation and evolu tion of proto-galaxies and the large scale s tructure in 



the early Universe (e.g., iKay et all 120021 ; IScannapieco et al.l . 120061 . and refer- 



ences therein). Although SN feedback is an important ingredient for galaxy 
formation, it occurs typically on length and time scales much smaller than grid 
spacings and time steps of numerical simulations. So such subgrid physics is 
often treated through phenomenological prescriptions in num erical simulations 
for g alaxy formation and large scale structure f ormation (ICen fc Ostrikerl . 



20061). Either 'the r mal' f eedback (e.g.. iKatzl . Il992l ) or 'kinetic' feedback (e.g., 
Navarro fc Whitd . Il993l ) with a feedback efficiency parameter < e < 1, 



which is the fraction of SN explosion energy deposited in the form thermal 
or kinetic energy, were tried with varying degrees of success. For a thermal 
feedback model, the radiative cooling must be turned off artificially for a pe- 
riod of time in order to prevent the de nse gas from re-radiating the deposited 
energy immediately. IKay et al.l (120021 ) showed that both thermal and kinetic 
feedback models can alleviate the so-called "catastrophic cooling" problem of 
galaxy formation. Thus in order to construct a physical subgrid prescription 
for SN feedback in numerical studies of galaxy formation it is necessary to 
estimate quantitatively the feedback efficiency parameters, e t h, tkim and e rad 
for thermal, kinetic and radiation energy, respectively. 



Massive stars often form in an associatio n, so multiple explosions of core- 



collap s e SNe can crea t e a la rge superb ubble (ITomisaka fc Ikeuchil . Il986l ; iTenorio et al. 



1990al ; ISmith fc Coxi . l200ll ). Multiple SN explosions in the proto-galactic en- 



vironment, in particular, may have played a significant role in the formation 
of halo globular clusters. In a self-enrichment model of metal-poor globular 
clusters, proto-globular cluster clouds are assumed to be enriched by multiple 
SN explosions from massive first-generation stars contained within the clouds 
and later seco nd-generation stars form in the cooled shell driven by the SN ex- 
plosions (e.g.. iParmetier et al.l . Il999l . and references therein). However, some 
have argued that the binding energy of a proto-globular cluster cloud may not 
be large enough to survive a few SNe, and that the homogeneity of metalicity 
within a globular clus ter indicates that pre -enrichment must have occurred 
within larger systems (jSearle fc Zinnl . Il978l ) . According to numerical studies 
of early structure formation, first objects might form at a redshift of ~ 30 with 
a typical mass of 10 6 M Q and first-generation stars enrich the protogalactic 
clouds before they merge into l arger clouds of 10 12 M via gravitational clus- 
tering (jBromm fc Larson! |2004| ). This favors the model that combines the pre- 
enrichment by first stars before the formation of protogalactic halos and the 
self-enrichment by next generation SNe within proto-globular cluster clouds. 
In this scenario the stability of the cooled shell and further star formation 
critically depends on how efficiently the SN energy is radiated away during 
the SNR evolution. 
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Thus in the present work we attempt to estimate feedback effects of core- 
collapse SNe, including multiple detonation, through numerical hydrodynamic 
simulations. We first calculate the creation of a bubble by stellar winds from 
multiple massive stars (N star = 1 — 50) during the precedent MS stage. The 
end results of the MS wind stage are adopted as the initial states of the SN 
explosion simulations. We then make quantitative estimates of the fractions of 
explosion energy transferred into the ISM in different forms, i.e., thermal, ki- 
netic, and radiation energy. Diffuse radiation emitted by the shock-heated gas, 
especially, ionizing radiation for hydrogen and helium atoms are calculated. 

In §2, the basic features of our numerical simulations, including radiative cool- 
ing and thermal conduction, are described. We then present the numerical 
results of wind bubble models in §3, and multiple supernova models in §4. A 
brief summary is given in §5. 



2 Numerical Method 



2.1 Hydrodynamic Equations 



We solve the hydrodynamic conservation equations of the ideal gas in one- 
dimensional spherical coordinate system, including radiative cooling and ther- 
mal conduction : 



0. 



d(pu) 1 d 2 2 

d(pe) 1 d 
dt r 2 dr 



dP 

dr ' 



r 2 (peu + Pu) 



(2) 
(3) 



where u is the radial velocity, e = P/\p(j — 1)] + u 2 /2 is the total energy of 
the gas per unit mass, A is the energy loss rate per unit volume, and S is the 
radial heat flux. A ratio of specific heats is assumed to be 7 = 5/3. The gas 
temperature is determined by T = fimnP/ipkB), where k^ is Boltzmann's 
constant, m# is the mass of hydrogen nucleus and p is the mean molecular 
weight. 



We use a grid-based, Eulerian hydrodynamics code in the spherical sym- 
metric geometry, based on the Total Variation Diminishing (TVD) scheme 
(IRyu et all Il993l ). The TVD scheme solves a hyperbolic system of gasdynam- 
ical conservation equations with a second-order accuracy. It has been widely 
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used for a variety of astrophysical problems (e.g., iBaek et all 120031 ). 



2.2 Radiative Cooling and Thermal Conduction 



We define the cooling rate as 
A(T,n H ) = L(T)n 2 H 



(4) 



where Uh is the number density of hydrogen nuclei. For the cooling rate coef- 
ficient, L(T), we adopt the non-equilibrium radiative cooling rate for an op- 
tically thin gas that is calculated by following the non-equilibrium collisional 
ioniz ation of the gas cooling from 10 8,5 K to 10 4 K under the isobaric condi- 
tion ([Sutherland fc Dopital . ll993l ). We consider cases with the solar metalicity 
( Z ) and with 1/10 of the solar metalicity (0.1 Z ). 

Although we calculate diffuse radiation emitted by the shock-heated gas, we 
do not explicitly follow the radiative transfer of diffuse radiation and pho- 
toionization of the gas. Instead we assume L(T) = for T < 10 4 K and set the 
minimum temperature at T min = 10 4 K. This is in effect equivalent to assum- 
ing thermal equilibrium at 10 4 K between radiative cooling and photoionization 
heating. 

After completing the hydrodynamic part, updating hydrodynamical quanti- 
ties, radiative co oling is app l ied to the thermal energy per unit volume, E, as 
a separate step (IBaek et all 120031 ) : 



E 



n+l 



E n ■ exp( 



At 



•'cool 



(5) 



where t coo i = E/A is the cooling time scale and At is the integration time 
step. This is designed to follow approximately the thermal history of the gas 
that cools radiatively and approaches the specified mini mum temperatur e in 
a time scale much shorter than the dynamical time scale (iLeVeque! . 119971 ). 



In some astrophysical plasmas such as solar corona and stellar winds ther- 
mal conduction can be important in heat transpor t. We in c lude the so-called 
classical Spitzer conduction for a fully ionized gas (jSpitzerl . 119621 ): 



S c i = -kV r T, 



1.84 x 10" 5 T 5 / 2 
In A, 



(6) 
(7) 
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where lnA c = 29.7 + Inn" 1 / 2 (T/10 6 K) and n e is the electron number density. 
This classical formula is based on the assumption that the mean free path of 
the electrons is short compared to temperature scale height. However, when 
temperature scale height is comparable or shorter than the electron mean 
free path, thermal conduction reaches the saturation limit and heat flux is no 
longer equal to /cV r T. In the saturation limit, heat flux is given approximately 
by 



S, 



sal 



-0.4( 



1 1 n e k B T, 



T\m r - 



where m e is the electron mass (jCowie fc McKed . 119771 ) . So we set S 
as local heat flux. 



mm 



sat] 



With the classical formula the thermal conduction term becomes a diffusion 
term, v r • (kv r T). This type of parti al differential equa tion can be solved 
implicitly by Cranck-Nicholson Scheme (IPress et all Il992l ). For the saturated 
formula, however, this cannot be applied. So we solve the thermal conduction 
term explicitly and the integration time step is chosen as 



At n = 0.1 -min(- 



E 



Ar 



(9) 



where c s = (7P//5) 1 / 2 is the sound speed. 

In grid-based Eulerian simulations a shock transition is typically resolved over 
several grid zones with a strong temperature gradient, although the physical 
shock thickness may be much smaller than the size of one grid zone. So un- 
physically rapid radiative cooling and heat transport can occur across the 
numerical shock transition zones. To avoid such numerical artifacts, we identify 
zones that undergo shocking process and turn off radiative cooling and heat 
transport by thermal conduction inside those zones. 

We neglect magnetic fields and dust in our calculations. If there exist dy- 
namically significant, regular magnetic fields in the ambient medium, the as- 
sumed spherical symmetry would brake down. Also coherent tangential mag- 
netic fields would substantially reduce the radial thermal conduction flux. 
Dust may play some important roles through radiation pressure on dust and 
dust cooling in the cooled bubble shell in a more realistic situation. Proper 
treatments of such processes require radiative transfer, chemistry and dust 
formation/destruction, which are beyond the scope of this study. However, if 
the bubble shell stays photoionized at 10 4 K by diffuse r adiation, as we assum e 



here, then the effects of dust are likely to be minimal. ISmith fc Coxi (120011 ). 
for example, included dust cooling in their simulations of multiple SNe and 
found that the dust cooling effects were only minor. 
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3 Development of Wind Bubbles 



3. 1 Wind Models 



We first consider the development of a wind bubble to study the pre-supernova 
circumstellar environment. Massive early-type stars emit a steady stellar wind 
with constant terminal velocity V w and mass-loss rate M w = dM w /dt. We as- 
sume that the progenitor star has an initial mass 15M and spends At w = 
4 x 10 6 yr in the main sequence phase. The MS wind is characterized with a 
mass-loss rate , M w = 2.5 x 10 -7 M yr -1 , and wind velocity, V w = 2000 km s _1 



([Smith! 120061 ). Here we consider only the MS wind, since the RSG wind pro- 
duces a small shell with the radius of a few pc inside the bubble created by the 
MS wind. The RSG wind shell would be swept up quickly by the SN ejecta 
during the early free expansion stage, so interactions of the ejecta with this 
shell can be ignored. 

We calculate the evolution of a bubble blown by collective stellar winds from 
multiple massive stars in an association. The luminosity of the collective wind 
is the key physical parameter that determines the properties and evolution of 
the bubble and it is given by 

L w = \m w V* ■ N star = N star • (3.17 x lO^ergs^ 1 )- (10) 



We select the number of stars, N sta r = 1, 5, 10, 30 and 50, because a typical OB 
association contains 10-100 early-type stars (see Table 1). In a realistic model, 
one needs to specify the initial mass function of the association to obtain a 
total wind luminosity. In our simplified model, however, L w , is parameterized 
with N s tar- So a larger value of L w can be regarded as a higher wind luminosity 
either from a more massive star or from multiple stars. 

We consider nap = 1 cm -3 as the fiducial density of the ambient ISM (WF 
models), where po = (2.34x 10~ 24 g)n.H;o. Two additional values are considered: 
n Hfl = 0.3 cm -3 (WL01 model) and 3 cm -3 (WH01 model). The ambient 
medium is assumed to be photoionized at T = 10 4 K. So its sound speed is 
c s = 14 km s _1 and the pressure is Pism = 2.3 x 10 _12 erg cm -3 . With WZ01 
and WZ10 models, we consider metal poor cases with 0.1 Z . 

The wind is realized numerically by injecting the gas flow with V w at the inner 
boundary (r w ) of the spherical simulation grid. The wind density at the inner 
boundary is given by p w = N star ■ [M w / (47rV^r^)]. We use N r = 3000 uniform 
spherical grid zones in 150 pc radius with a spatial grid spacing of Ar = 0.05 
pc. Numerical convergence is often an important issue for hydrodynamic sim- 
ulations with radiative cooling, because the cooling rate has a strong density 
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dependence (i.e., A oc p 2 ). So in order to test numerical convergence of the 
simulations WF01 model is calculated with several different grid spacings. Fig. 
1 shows the structure of the bubble in the simulations with 3 different grid 
spacings. The contact discontinuity (CD) between the shocked wind and the 
shocked ISM is around 23 pc, while the location of the conduction front (CF) 
is different in the three simulations. The size of the hot cavity (~ 20pc) is 
smaller than the so-called Field length A f = \kT/ (njrL)] 1 ^ 2 & 56 p c, where 
T ~ 10 6 K and n H ~ 0.01 cm" 3 are used flBegelman fc McKed . Il990l ). So heat 
transport by thermal conduction is dominant over radiative cooling inside the 
cavity. Heat transport from the hot shocked wind gas to the warm shocked ISM 
reduces the temperature of the shocked wind gas, triggering radiative cooling 
and compression. These processes are not followed accurately enough to be 
converged in the simulations, so the density and temperature distributions 
are significantly different inside the cavity among these simulations. However, 
the pressure distribution is very similar in the three calculations and so the 
dynamical evolution of the wind bubble shell seems to be converged. Radial 
positions of the CD and outer compression wave seem to be converged within 
a few percent, while the total energy inside the simulation volume seems to 
be converged within ~ 10 %. This fast convergence is achieved, because the 
density compression factor at the shell is not very large in the simulations 
{p shell / Po ^ 10). We also find that the differences in the density and temper- 
ature structures inside the cavity among the simulations with different grid 
spacings do not affect significantly the evolution of SN explosions that will be 
described in 54.2. 



3.2 Evolution of Wind Bubbles 



A typical bubble has four distinct regions from the central star outward: 1) 
freely expanding wind with the density profile, p oc r~ 2 , 2) a wind termination 
shock facing inward and the shocked wind region, 3) a contact discontinuity 
(CD) between the shocked wind region and the shocked ISM regio n, 4) a 



forward shock propagating into the ambient ISM ([Weaver et al.l . 119771 ). While 
the termination shock is strong, a compression wave may form, instead of 
the outer forward shock, depending on the ratio of wind ram pressure to the 
ambient pressure. 

Fig. 2 shows the gas density and temperature profiles at t = 10 6 years in 
the models listed in Table 1. Although the outermost structure is initially 
developed as a weak shock for all models, it becomes a compression wave with 
the propagation speed, v < c s in some models, as the bubble expands out into 
the ambient medium. For WF01 model, for example, the pressure increases 
only by a factor of 1.5 across the compression wave. Even for WF50 model 
the outer shock speeds is only v s ~ 20 km s _1 . Also the shock transition 
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is much wider than that typically found in adiabatic simulations using the 
same TVD code. In WF50 model the inner termination shock is at 10 pc, the 
conduction front at ~ 25 pc, the CD at 41 pc, and the outer shock wave at 
Rb s ~ 45 pc. Without thermal conduction, the hot cavity inside the CD would 
be almost isothermal with T « ]_g".6-7.8 j£ g u ^- h ea t transport by thermal 
conduction enhances radiative cooling near the conduction front, leading to a 
strong temperature gradient inside the cavity. 

Near the CD, the gas density changes by more than two orders of magnitude, 
so the gas cools quickly down to T min = 10 4 K and forms a shell. The thickness 
of the shell is rather thick (AR bs ~ 0.3i?f, s ). The continuous injection of wind 
energy seems to be radiated away efficiently near the conduction front. For a 
spherical shell with inner radius, Rb s ,i = aRbs, outer radius, Rb s , and uniform 
density, the shell density can be derived from the condition that the shell mass 
is equal to the swept-up mass: 

^^(1-a 3 )- 1 . (11) 



So the density enhancement factor is small for thick shells found in our simu- 
lations. For WF01 model, for example, a ~ 0.7 and n sheU /n H fi « 1.5. 

From Fig. 2 we c an see that th e size of bubble increases with N star and de- 



creases with n H) Q. I Weaver et al.l (119771 ) showed that the bubble would evolve 



self-similarly and the outer radius of the shell would grow as R bs oc £j/ 5 Po 1//5 t 3 / 5 , 
if radiative cooling and the ambient pressure were negligible. They also derived 
a slightly modified time dependence of Rb s oc t 58 by numerical integrations 
of simplified conservation equations where an approximate cooling rate and 
an ambient pressure with Tism =8000K were adopted. In our simulations the 
ambient pressure is not negligible and the outward moving front is not a shock, 
so the self-similarity seems to be broken. 

Fig. 3 shows the evolution of Rb s and swept-up shell mass, M bs = (47r/3)Rf s p , 
for all the models in Table 1. For our models with L w = N star (3.17x 10 35 ergs _1 ) 
and nH,o = 1 cm -3 , the radius and mass of wind-blown bubble shell can be 
approximated by 

R bs « (22 P c)iVr (y^) - 4 , (12) 



and 



Mbs « (1.6 x 10 3 M )C(^) L2 (13) 
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These relations are slightly different from what IWeaver et al.l (119771 ) derived 
through simplified numerical integrations described above. 



Right panels of Fig. 3 show how Rb s and depend on uh,o or the metalicity 
Z. These quantities do not scale with the ambient density in a simple way, 
because the radiative energy loss depends not only on the gas density but 
also on the detailed history of how heat is transported from the hot cavity to 
the warm shell via thermal conduction. It seems that for nu,o < 1 cm -3 the 

— 1/3 

swept-up mass is nearly independent of the ISM density, while Rb s oc n H q . For 
high density medium (njj,o ^ 1 cm -3 ), the cooling time scale becomes much 
shorter than the dynamical time, and so the numerical simulation becomes 
increasingly difficult. For WH01 model with Uh,o = 3 cm -3 , M& s is about 
3 times that of WF01 models. However, this trend cannot be generalized 
to higher density cases because of nonlinear behavior of heat transport and 
radiative cooling. For WZ01 with lower metalicity (0.1 Z ) the cooling is less 
efficient and so the results are similar to those of WL01 model with lower 
ambient density. 

Bottom panels of Fig. 3 show the fraction of energy lost due to radiative cool- 
ing. The radiation energy is calculated as the energy lost from the simulation 
volume, E rad (t) = L w ■ t + E tot ,i - E tot (t), where E tot = E k + E th is the total 
energy inside the simulation volume and E tot ^ is its initial value. Here the 
thermal energy E t h excludes the initial thermal energy of the ambient gas. 
The energy loss due to radiative cooling is significant: about 60 % for WL01 
and WLZ01 models and about 60-80 % for WH01 and WF01-WF50 models. 
For WF models the fraction of E rad is reduced somewhat for higher L w (larger 
N s tar), because higher power wind creates hotter cavities. 



4 Multiple Supernova Explosion 



4-1 Multiple Supernova Models 



Next we calculate the evolution of multiple supernovae exploded inside the pre- 
existing wind bubble created by the MS winds. So we take the final structures 
of the wind bubble simulations at t — 4 x 10 6 yrs, described in the previous 
section, as the states of the circumstellar medium prior to SN explosions. 

Each SN explosion is characterized by the explosion energy, E e j = 10 51 erg, 
and the ejecta mass, M e j = 10 M . Then the ejecta can be approximated 
by a uniform density core (pi) and an outer envelope with a steep power-law 
density profile {pi e x r~ n ) , expanding freely with the radial velocity, v = r/t 



( Chevalier fc Liang . 1989 ). For a multiple explosion where N S n supernovae 
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detonate simultaneously, i.e., total explosion energy E$n = N$NE e j, the den- 
sity profile of the ejecta is defined by 



Pi — Ft for v < v t , 



P2 = Ft [—) for v > v t, 



(14) 
(15) 



where 



N SN [3(n - 3)M ej ] 5/2 
47rn [io(n-5)E ej ] 3/2 ' 



and 



v t = 



10 (n - 5) E ej 
T (n - 3) M ej 



1/2 



(17) 



We adopt n = 10, u t = 3.16 x 10 3 km s" 1 and F = N$n~ (8.14 x 10 7 g cm' 3 s 3 ). 
We note that overall results of our calculations do not depend sensitively on 
the specific values of M ej and E e j adopted. As shown below, the evolution of 
a superb ubble is controlled mainly by the total explosion energy Esn rather 
than E e j. Also a SNR should lose its memory of ejecta mass after the swept-up 
mass becomes much greater than M e j. 

Table 2 shows the multiple SN explosion models considered here. The number 
of SNe is given by N$ni and varies from 1 to N star , since N star is the total 
number of possible SNe inside a given bubble. As in the wind bubble stage, 
we consider three values of the ISM density: nn,o = 1 cm -3 as a fiducial value 
(SNF models), 0.3 cm" 3 (SNL01) and 3 cm" 3 (SNH01). Since SN explosions 
can occur sequentially with a time interval, we also consider SNF02a model in 
which a second SN detonates 10 5 years after the first SN. This will be compared 
with SNF02b model in which two SNe explode simultaneously. Again SNZ01 
and SNZ10 models have a lower metalicity, 0.1 Z Q . 

We use the same grid setup as the wind simulations: i.e., N r = 3000 uniform 
grid zones in 150 pc radius. Again two additional simulations are carried out 
for SNF01 model with smaller grid spacings (N r = 6000 and 12000) to test 
numerical convergence. Radial position of the outer shock is converged within 
0.5 %, while the total energy is converged within 5 %. The simulations are 
terminated at tf — 10 6 years. 
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4-2 Evolution of Supernovae Remnants 

Fig. 4 shows the radial profile of the remnant for the single SN case (SNF01 
model). Initially, the SN profile described by Eqs. (14)- (17) is superimposed 
onto the final output of the corresponding wind bubble model (WF01): the 
wind termination shock is at 3 pc and the hot cavity extends out to 25 pc (solid 
lines in Fig. 4). The compression wave has weakened and the shell thickness is 
large (a ~ 0.43), so the density of swept up shell, n s h e ii = 1.1 cm -3 , is slightly 
higher than the ambient density (see Eq. [11]). The SN blast shock wave 
passes through the pressure barrier produced by the wind termination shock 
and travels quickly through the hot cavity. The forward shock first interacts 
with the bubble shell at ~ 10 4 years, pushing out the inner shell boundary and 
generating two secondary waves, i.e., a transmitted shock propagating outward 
into the bubble shell and a weak reflected compression wave moving inward 
into the ejecta. A reverse shock forms in the SN ejecta around 5 x 10 3 years 
and is reflected at the center around 7.5 x 10 3 years. Afterward it becomes a 
traveling compression wave, because the shock heated gas inside the cavity is 
hot. When this wave runs into the shell at ~ 10 5 years, the shell is heated for 
the second time. 

By the termination time (10 6 years) the shell is pushed out by the SNR to 
R ss =67 pc with inner radius R ss ,i =54 pc and the swept up shell gas has 
T = 10 4 K and n s h e u ~ 4 cm -3 , not much higher than the ambient medium. 
On the other hand, the cavity is filled with the hot gas of T 10 6 K and 
uh ~ 0.01 cm~ 3 . The outer shock weakens with expansion and slows down to 
~ c s . 

Fig. 5 compares the two different models of double explosions, i.e., sequential 
(SNF02a) and simultaneous (SNF02b) explosion models. In SNF02a model, 
the second SN explosion is loaded 10 5 years after the first SN explosion. Three 
snapshots are shown: 1) At 6.0 x 10 4 years the outer shock in SNF02b model 
has expanded further out than the shock in SNF02a model. 2) At 1.2 x 10 5 
years, a little after loading of the second SN in SNF02a model, the second SN 
shock expands out to 35 pc, still well inside the shell at 42-46 pc generated by 
the first SN. 3) At 2.0 x 10 5 years the SN-driven shell has expanded slightly 
further out (~50 pc) in SNF02b model than in SNF02a model. However, the 
difference is small enough, so we say the bubble size is similar for the two 
models. 

The left panels of Fig. 6 also show the comparison of the two double explo- 
sion models. The outer shock radius and shell mass of the two models are 
very similar. Here E t h and are total thermal and kinetic energy inside the 
simulation volume, respectively, excluding the initial thermal energy of the 
ambient gas. The time evolution of E th and E k shows that the initial kinetic 
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energy of the SN ejecta is transferred to the thermal energy gradually in the 
beginning as the blast wave expands, but rather abruptly at ~ 10 4 years when 
the blast wave runs into the preexisting wind-blown shell. The thermalization 
process happens again at later times in SNF02a model, when the second SN 
detonates at 10 5 years and when its blast wave arrives at the first SN shell at 
~ 10 5 ' 2 years. Interestingly, the latter event occurs at the similar time when 
the aforementioned compression wave hits the shell in SNF02b model, causing 
the second period of active radiative cooling. 

Most of thermal energy is lost through radiative cooling. Here the radiation 
energy is calculated by E ra( i{t) = E to t,i — E tot (t). We note that the initial total 
energy, E to t,i, includes both the SN energy and the energy of the preexisting 
wind bubble, so (E tot +E rad )/ E SN is greater than one. At the termination time 
of our simulations, the final energy budget in the two models are similar: 6% in 
thermal energy, 7% in kinetic energy, and 90% lost through radiation energy. 
Figs. 5-6 imply that, for given total explosion energy, the different fashion of 
multiple explosion, i.e., either sequential or simultaneous detonation, does not 
make significant differences in the size and mass of the final bubble shell. 

The right panels of Fig. 6 compare the models with different ambient density 
and different metalicity. At lower density or with lower metal abundance, the 
bubble size is bigger and E rad is smaller, compared to higher density models, 
but the shell mass is similar for all the models considered here. 

Simultaneous explosions with 1, 3, 5, and 10 SNe inside a wind bubble gener- 
ated by 10 massive stars are considered as SNFlOa-d models. Initial conditions 
are taken from WF10 model simulation: the cavity is bounded by the CD at 56 
pc and the outer compression wave is at 82 pc. The upper two panels of Fig. 7 
show the bubble structure at 10 6 years for SNFlOa-d models. More energetic 
explosion pushes the shell further out and produces a faster shock, leading to 
a larger and denser shell. The lower panels of Fig. 7 show the profiles of SNRs 
in a uniform ISM (nn,o — 1 cm~ 3 ) for comparison. It shows that SN-driven 
bubbles are larger and hotter, if SNe explode inside a preexisting bubble. For 
example, the size of hot bubble for SNFlOa-d models is 33-55 % larger than 
that for SNe inside a uniform ISM, resulting in 2.4 - 3.7 times bigger volume 
of hot gas. Also the cavity temperature is 2-3 times higher in SNFlOa-d mod- 
els. This implies that the volume filling factor and the temperature of the hot 
bubbles generated by core-collapse SNe would be underestimated by a factor 
of 2-4, if the preexisting bubble is not properly accounted for. 

The left panels of Fig. 8 show the evolution of the bubble properties for 
SNFlOa-d models, while the right panels compare the SNZlOa-d models with 
lower metalicity. The radial position of the outermost shock increases as R ss oc 
^0.6-0.8 kgforg the arrival of the blast wave at the shell, and afterward it in- 
creases very slowly, R ss oc t - 1 " 015 . The size of SN-driven bubbles depends on 
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both the SN explosion energy and the size of pre-existing wind bubbles. Ac- 
cording to the numerical results shown in Fig. 8, R ss and M ss depends rather 
weakly on E SN} i.e., R ss w (85pc)JV[$ and M ss « (10 4 - 8 M )iVg#. 



Again the evolution of E t h, Ek, and E ra d indicates the complex thermalization 
history via shocks. The first major episode of radiative cooling occurs when 
the SN blast wave first hits the wind bubble at t ~ 10 4 ' 5 years. When the 
compression wave generated by the reverse shock runs into the shell at t ~ 10 5,4 
years, the shell is heated and compressed, leading to the second episode of 
rapid cooling. In SNFlOd model this behavior is most distinct. For SNFlOa-d 
models, about 80-90 % of E S n is lost via radiative cooling within 10 6 years, 
while only 10 % remains as thermal energy and <10 % as kinetic energy. 
For SNZlOa-d models with lower metal abundance the radiative cooling is less 
efficient compared to SNFlOa-d models, so about 70 % is lost, 20-40 % remains 
as thermal energy and 10 % as kinetic energy. Since the shell has a significant 
kinetic energy at that time, it would continue to expand slowly. Exploration of 
later evolutionary stage requires more realistic treatments of radiative transfer, 
photoionization/heating, and non-equilibrium radiative cooling below 10 4 K, 
which is beyond the scope of this study. 



In the self-enrichment model of globular cluster formation, the shell of gas 
swept up by the blast wave from multiple SNe must cool and be confined 
gravitationally within the pare nt proto-globula r cloud in order for the e nrich- 



ment process to be effective (IMorgan fc Lakd . Il989l ; iParmetier et all 11999 



Parmentierl . 120041 ) . Thus the success of this model lies on how efficiently the SN 



3 and 
« 50 



explosion energy is radiated away to form gravitationally bound system. Ac- 
cording to Fig. 8 (SNZ10 models), in an average ISM with uh,q = 1 cm" 
Pism = 2.3 x 10 _12 erg cm -3 , massive stars in an association with N star 
would generate a bubble shell with R ss ps 120 pc and M ss ps 10 5 5 M . In- 
cidentally this mass is similar to the characteristic mass of globular clusters. 
If 10 % of explosion energy is left as the kinetic energy of the shell, then 
Ek ~ 5 x 10 51 ergs. With the gas temperature of T = 10 4 K, the total ther- 
mal energy of the photoionized shell is E t h ~ 10 51 ergs. But the gravitational 
energy, \Eq\ ~ GM^ S / R ss « 10 50 ergs, is too small to confine the shell gravita- 
tionally. Thus the self-enrichment model would not work inside such low den- 
sity clouds, but may operate in much higher density environment. If we blindly 
take the simple scaling relations shown in Figs. 6 and 8, i.e., R ss oc njj 1 ^ , but 
M ss p^ constant, then the gravitational energy of the shell would increase 

1/3 

with the ambient density as \Eg\ oc n# . In a much denser cloud higher am- 
bient pressure would also help the shell to be confined, possibly resulting in a 
gravitationally bound system. 
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4-3 Diffuse Radiation Emitted by Supernova Driven Bubbles 



As shown in the previous section, most of mechanical SN energy is transformed 
into thermal energy of the bubble shell and then lost via diffuse radiation. So 
radiative feedback is a significant part of SN feedback into the surrounding 
ISM. In order to study the photoionization/heating feedback of SN explosions 
we calculate the diffuse radiation emitted by the bubble. 

First, we generate both continuum and line emissivities, f v = 4nj u (T) /n 2 H 
(ergs cm 3 s^eV" 1 ) and f Hne = Anj line (T)/n 2 H (ergs cm 3 s" 1 ), for n H = 1 cm" 3 
as a function of gas temperature for photon energy E p h(eV) —hv— [E min , E max ] = 
[1.36 x 10~ 4 , 10 8 ], using CLOUDY90 code (Ferland et al. 1998). Collisional ion- 
ization equilibrium is assumed and no background radiation field is adopted. 
The emissivity tables are generated for two values of metalicity, Z Q and 
0.1 Z . Volume integrated spectrum of the diffuse radiation emitted by the 
bubble during the integration time t is calculated by 

n 2 H (T,r)dr, (18) 



Fu{t) 



dt 



U{T) + 



fi 



line 



AE, 



ph 



where F v is given in units of ergs eV _1 and line emissivities divided by the 
width of a given photon energy bin are added to the corresponding continuum 
emissivity. In practice, only V(T,t) = An j dt J n 2 H (T)r 2 dr is calculated during 
simulations and then F u (t) = V(T,t) ■ [f v (T) + fu ne / AE ph ] is calculated in a 
post-processing step. 

The upper panel of Fig. 9 shows the evolution of F v (t) for a single SN model. 
Most of hard X-ray photons with E ph > 1 keV are emitted early before the 
blast wave runs into the wind bubble shell (t < 8 x 10 3 years). Afterward, 
the fluxes of lower energy photons increase with time, as the SN energy is 
deposited as heat into the bubble shell and radiated away mainly in UV and 
optical ranges. The lower panel of Fig. 9 compares F v for SNF01 and SNZ01 
models integrated to the final simulation time. In SNZ01 model metal lines are 
weaker, but H and He recombination continuum radiation is stronger, because 
the temperature is higher due to lower cooling rate. 

We also calculate the integrated energies and photon numbers of diffuse radi- 
ation that can ionize hydrogen and helium atoms as follows: 



£Ei = J dE ph F v (tf), 



(19) 
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E^ rn a 

I 



dE, 



ph' 



Mtf) 



E, 



ph 



(20) 



where Ei = 13.6, 24.6, and 54.4 eV is the ionization energy for H, He, and 
He + , respectively. These quantities are summarized in Table 3. For SNF01 
model, the total number of H ionizing photons is $i 3 . 6 w 10 61 and the total 
energy carried away with those photons is £13.6 ~ 0.55^^. The total number 
of H ionizing photons scales roughly with the SN explosion energy as $13.6 ~ 
10 61 AW 

Diffuse radiation is enhanced overall in the models where a denser and hotter 
shell is produced. Comparing the results of SNL01, SNF01 and SNH01 models 
in Table 3, we can see that more diffuse photons are produced at denser 
environments. A SN with higher explosion energy generates a faster blast 
wave and produces a denser and hotter shell, emitting harder diffuse radiation. 
The results of SNFlOa-d models show such trends. On the other hand, the 
comparison of SNF02a and SNF02b models indicate that the simultaneous 
double explosion generates more ionizing photons than the sequential double 
explosion, although the total radiation energy (E ra d) is similar in the two 
models. SNZ models with lower metalicity produce less ionizing radiation and 
more non-ionizing radiation, compared to the models with solar metalicity. 

Fig. 10 shows the evolution of for SNFlOa-d and SNZlOa-d models. Most 
of ionizing photons are emitted during 10 4 5 — 10 5 years when the SN shock 
propagates into the bubble shell. For models with N S n = 10 ionizing photons 
are emitted actively again when the reflected compression wave runs into the 
shell around 10 53 years (see also Fig. 8). 



5 Summary 



We study multiple supernova (SN) explosions inside a preexisting cavity blown 
up by stellar winds from massive progenitor stars through numerical hydro- 
dynamic simulations. Calculations are performed with a grid-based Eulerian 
code in one dimensional spherical symmetry, including radiative cooling and 
thermal conduction. Heat transport by thermal conduction is important in- 
side the hot cavity, while radiative cooling is dominant in the bubble shell. 
We consider models with a wide range of physical parameters: the ambient 
density, 0.3 < nn,o < 3 cm -3 , the metalicity, Z Q and 0.1 Z Q , the number of 
massive stars, N star = 1 — 50, and the number of supernovae, N$n = 1 — 10. 
In order to study the radiative feedback effects, we also calculate the ionizing 
diffuse radiation emitted by the SN-driven bubble. 



16 



We first calculate the development of a bubble driven by collective winds from 
N star massive stars with L w = N star ■ (3.17 x lO^ergs" 1 ). Characteristic struc- 
tures of wind bubbles are developed initially: an inward facing termination 
shock, a conduction front near the contact discontinuity, and a forward ex- 
panding shock. The outermost shock is not strong and becomes a compression 
wave, because the wind ram pressure decreases as the bubble expands. At the 
end of the main sequence (~4x 10 6 years), the outer boundary of the bubble 
shell (i.e., radius of the compression wave) is R bs ~ (60pc)iV° t ^, and the inner 
boundary is R bs ,i ~ (0.7 — 0.9)R bs for % i0 = 1 cm -3 . So the shell is rather 
thick and the enhancement factor of the shell density is small, i.e., n s h e ii/nH,o — 
1.5 - 3.7. The mass of the bubble shell is M bs « (2 x 10 4 M )A^ 4 . 

The evolution of a supernova remnant inside a wind bubble is much more 
complex than that in a uniform medium. It depends on the structure and size 
of the preexisting bubble as well as the total explosion energy, E S n = A^-IO 51 
ergs. The forward shock travels quickly through the cavity and runs into the 
bubble shell, resulting in a transmitted forward shock. On the other hand, the 
reverse shock in the ejecta is reflected at the center and becomes a compression 
wave that travels outward. When this wave runs into the shell, the shell is 
heated for the second time, followed by the second period of active radiative 
cooling. The wave bounces off the shell and moves inward with a much reduced 
amplitude. 

At the termination time of the simulations (tf — 10 6 years), the outer shock 
radius is R ss ~ 67 pc for a single SN inside a bubble blown by one massive 
star (SNF01 model), while R ss ~ 85 pc for a single SN inside a bubble blown 
by 10 massive stars (SNFlOa model). In comparison, the shock radius for a 
single SN in a uniform ISM of the same density at the same time is about 
60 pc. The radii of the hot cavity and the bubble shell can be bigger up to 
50 % and the cavity temperature is higher by 2-3 times for SN explosions 
in a wind bubble considered here, compared to the SN explosions inside a 
uniform ISM. For multiple SN explosions inside a bubble blown by 10 massive 
stars, the size and mass of the shell increase rather weakly with the total 
explosion energy as R ss (85pc)A^ and a mass of M ss (10 4 - 8 M )iV^. 
For N S n = 50, for example, the SN-driven shell would have typically R ss ~ 125 
pc and M ss ~ 10 5 - 3 M for the average ISM. 

Kinetic energy of the SN is transfered to thermal energy via various shocks. 
About 20-60 % of the explosion energy is radiated away after the SN shock 
hits the wind bubble shell, and then 20-40 % is lost after the compression 
wave crushes into the shell later. Although the final energy budget at the 
termination time of our simulations depends on both the wind bubble structure 
and Esn, in most cases with solar metalicity, 2±w/ ' E$n ~ 0.8—0.9, E t h/EsN ~ 
0.1 and Ek/EsN 0.1. Ionizing photons of diffuse radiation account for up 
to £ i3.6/ 'Esn ~ 0.55, while the total number of ionizing photons scales as 
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$13.6 ~ i0 61 N S N- For the cases with 0.1 Z , the radiative loss is reduced a little 
bit, i.e., E ra d/EsN ~ 0.7, but the fraction of ionizing radiation is reduced a 
factor of 2-3, compared to the models with solar metalicity. So a larger fraction 
of energy is lost via non-ionizing radiation at metal-poor environments. 

We suggest physically correct prescriptions for SN feedback in numerical sim- 
ulations of galaxy formation should reflect the natures described in this study. 
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Table 1. Wind Model Parameters 



model 


nH,o (cm -3 ) 


star 


Z/Z Q 


WF01 


1 


l 




WF02 


1 


2 


1 


WF05 


1 


5 


1 


WF10 


1 


10 


1 


Wr oU 


i 
l 


on 

oU 




WF50 


1 


50 




WL01 


0.3 


1 




WH01 


3 


1 




WZ01 


1 


1 


0.1 


WZ10 


1 


10 


0.1 



Table 2. Multiple Supernova Explosion Models 



model 


riH,o (cm -3 ) 


N s tar 


Nsni 


NsN2 


Z/Z Q 


SNF01 


1 


1 


1 





1 


SNF02a 


1 


2 


1 


1 


1 


SNF02b 


1 


2 


2 





1 


SNFlOa 


1 


10 


1 





1 


SNFlOb 


1 


10 


3 





1 


SNFlOc 


1 


10 


5 





1 


SNFlOd 


1 


10 


10 





1 


SNL01 


0.3 


1 


1 





1 


SNH01 


3 


1 


1 





1 


SNZ01 


1 


1 


1 





0.1 


SNZlOa 


1 


10 


1 





0.1 


SNZlOb 


1 


10 


3 





0.1 


SNZlOc 


1 


10 


5 





0.1 


SNZlOd 


1 


10 


10 





0.1 
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Table 3. Ionizing Photon Number and Energy of Diffuse Radiation 



model log($i 3 . 6 ) log($24.e) log($ 5 4.4) £13.6/ 'E S N £24.6/ E S N £54.4/ E S N 



SNF01 


61.2 


60.4 


59.8 


5.5(-l) 


1.6(-1) 


7.7(-2) 


SNF02a 


61.3 


60.1 


59.3 


2.8(-l) 


4.1 (-2) 


1.4(-2) 


SNF02b 


61.5 


60.7 


60.1 


5.5(-l) 


1.9(-1) 


7.9(-2) 


SNFlOa 


60.7 


59.8 


59.2 


1.6(-1) 


4.5(-2) 


2.3(-2) 


SNFlOb 


61.5 


60.6 


59.6 


3.5(-l) 


7.4(-2) 


1.9(-2) 


SNFlOc 


61.8 


61.0 


60.1 


4.6(-l) 


1.2(-1) 


3.6(-2) 


SNFlOd 


62.2 


61.5 


60.9 


5.8(-l) 


2.2(-l) 


9.9(-2) 


SNL01 


61.0 


60.2 


59.6 


3.6(-l) 


1.2(-1) 


4.7(-2) 


SNH01 


61.4 


60.6 


59.9 


8.6(-l) 


2.9(-l) 


l.l(-l) 


SNZ01 


60.8 


60.1 


59.2 


2.2(-l) 


6.8(-2) 


1.9(-2) 


SNZlOa 


60.3 


59.3 


58.6 


6.4(-2) 


1.3(-2) 


5.3(-3) 


SNZlOb 


61.2 


60.5 


59.3 


1.7(-1) 


5.4(-2) 


8.8(-3) 


SNZlOc 


61.5 


60.8 


59.9 


2.4(-l) 


7.2(-2) 


1.9(-2) 


SNZlOd 


61.8 


61.1 


60.3 


2.3(-l) 


7.4(-2) 


2.5(-2) 
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Fig. 1. Gas density, temperature, and pressure profiles at t = 2.5 x 10 6 years 
in WF01 model calculated with N r = 3000 (solid lines), 6000 (dotted), and 12000 
(dashed) radial grid zones. The density is given in units of p a = 2.34 x 10 _24 g cm -3 , 
and the pressure in units of P Q = 9.36 x 10 _8 erg cm~ 3 . 
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Fig. 2. Upper panels: Gas density and temperature profiles at t = 10 6 years in WF01 
(solid lines), WF10 (dotted), WF30 (dashed), and WF50 (long dashed) models. 
Lower panels: Gas density and temperature profile at t = 10 6 yr for WF01 (solid 
line), WL01 (dotted), WH01 (dashed), and WZ01 (long dashed) models. The density 
is given in units of p a = 2.34 x 10~ gem -3 . 
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Fig. 3. Left panels: Radial position of the outer compression wave (top), swept-up 
mass (middle) of the wind-blown bubble shell, and the faction of energy lost due 
to radiative cooling for WF01 (solid lines), WF03 (dotted), WF05 (dashed), WF10 
(long dashed), WF30 (dot-dashed), and WF50 (dot-long dashed) models. Right 
panels: Same as the left panels but for WF01 (solid lines), WL01 (dotted), WH01 
(dashed), and WZ01 (long dashed)models. The radiation energy E ra( i is given in 
units of the total energy injected by winds. 
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Fig. 4. Time evolution of the supernova remnant from a single SN (SNF01 
model) at t = (solid lines), 10 4 years (dotted), 10 5 years (dashed), 5. x 10 5 
years (long dashed), and 10 6 years (dot-dashed). The normalization constants are : 
p = 2.34 x 10" 24 gcm- 3 , u Q = 2000 km s _1 , and P Q = 9.37 x 10^ 8 erg Cm - 3 . 
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Fig. 5. Time evolution of the supernova remnant for double SN explosion cases. 
Solid lines are for the sequential explosion model (SNF02a), while dashed lines are 
for the simultaneous explosion model (SNF02b). Three profiles correspond to the 
remnants at t = 6 x 10 4 yr, 1.2 x 10 5 yr, and 2 x 10 5 yr, from inside to outside, 
respectively. The normalization constants are the same as Fig. 3. 
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Fig. 6. Left panels: Radial position of the outer shock, swept-up mass of the SN- 
driven shell, thermal energy, kinetic energy, and radiation energy (from top to bot- 
tom, respectively) are shown for SNF02a (solid lines) and SNF02b (dashed lines) 
models. Right panels: Same as the left panels except that the results for are shown 
for SNF01 (solid lines), SNL01 (dotted lines), and SNH01 (dashed), and SNZ01 
(long dashed) models. 



27 




Fig. 7. Upper panels: Supernova explosions inside a bubble created by 10 massive 
stars (WF10 model). SNFlOa (solid lines), SNFlOb (dotted), SNFlOc (dashed), 
and SNFlOd (long dashed) are shown at t = 10 6 years. Lower panels: Super- 
nova explosions inside a uniform ISM of njjfi = 1 cm -3 . The explosion energy 
is Esn = NsnW 51 ergs, where N$n = 1 (solid lines), 3 (dotted), 5 (dashed), and 
10 (long dashed). Here p = 2.34 x 10 _24 gcm~ 3 . 
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Fig. 8. Left panels: Radial position of the outer shock wave, swept-up mass of the 
SN-driven shell, thermal energy, kinetic energy, and the energy lost due to radiative 
cooling (from top to bottom, respectively) for SNFlOa (solid lines), SNFlOb (dot- 
ted), SNFlOc (dashed), and SNFlOd (long dashed) are shown. Right panels: Same 
as the left panels except that Z = O.IZq is assumed. 
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Fig. 9. Upper panel: Integrated spectrum of diffuse radiation emitted by the bubble 
for SNF01 model at t = 10 years (red), 5 x 10 4 years (blue), and 10 6 years (green). 
Lower panel: Integrated spectrum of diffuse radiation emitted by the bubble in 
SNF01 model with Z = Z Q (red) and SNZ01 with Z = 0.1Z Q (blue) models at 
t = 10 6 years. 
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Fig. 10. Time evolution of integrated ionizing energies of diffuse radiation emitted 
by the bubble in SNF10 models with Z = Z & (left panels) and SNZ10 models 
with Z = O.IZq (right panels). The explosion energy is E$n = -^stvIO 51 ergs, 
where N$n = 1 (solid lines), 3 (dotted), 5 (dashed), and 10 (long dashed). Here 
p = 2.34 x 10 _24 gcm~ 3 . 
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